The Chirp z-Transform Algorithm
and Its Application
S By LAWRENCE R. RABINER, RONALD W. SCHAFER,
and CHARLES M. RADER*
(Manuscript received November 21, 1968)
We discuss a computational algorithm for numerically evaluating the z-
transform of a sequence of N samples. This algorithm has been named the
chirp z-transform algorithm. Using this algorithm one can efficiently
evaluate the z-transform at M points in the z-plane which lie on circular
or spiral contours beginning at any arbitrary point in the z-plane. The
angular spacing of the points is an arbitrary constant; M and N are
arbitrary integers.
The algorithm is based on the fact that the values of the z-transform on a
circular or spiral contour can be expressed as a discrete convolution. Thus
one can use well-known high-speed convolution techniques to evaluate the
transform efficiently. For M and N moderately large, the computation
time is roughly proportional to (N + M) log 2 (N + M) as opposed to
being proportional to N-M for direct evaluation of the z-transform at M
points.
Applications discussed include: enhancement of poles in spectral analysis,
high resolution narrow-band frequency analysis, interpolation of band-
limited waveforms, and the conversion of a base 2 fast Fourier transform
program into an arbitrary radix fast Fourier transform program.
I. INTRODUCTION
In dealing with sampled data the 2-transform plays the role which
is played by the Laplace transform in continuous time systems. One
example of its application is spectrum analysis. The computation of
sampled ^-transforms, which has been greatly facilitated by the fast
Fourier transform algorithm, is further facilitated by the "chirp
z-transform" algorithm described in this paper. 1 - 2
* Mr. Rader is with Lincoln Laboratory, Massachusetts Institute of Tech-
nology, Lexington, Massachusetts. Lincoln Laboratory is operated with support
from the U. S. Air Force.
1249
1250 THE BELL SYSTEM TECHNICAL JOURNAL, MAY-JUNE 1969
The z-transform of a sequence of numbers x n is denned as
X(z) = 2 x n z~ n , (1)
»- — CO
a function of the complex variable z. In general both x n and X(z)
could be complex. It is assumed that the sum on the right side of equa-
tion (1) converges for at least some values of z. We restrict ourselves
to the z-transform of sequences with only a finite number N of non-
zero points. Therefore, we can rewrite equation (1) without loss of
generality as
AT-l
X(z) = Z *** (2)
i =
where the sum in equation (2) converges for all z except z = 0.
Equations (1) and (2) are similar to the defining expressions for
the Laplace transform of a train of equally spaced impulses of magni-
tudes x„. Let the spacing of the impulses be T and let the train of
impulses be
Y,x n b{t-nT).
Then the Laplace transform is
n
which is the same as X (z) if we let
z = e' T . (3)
For sampled waveforms the relation between the original waveform
and the train of impulses is well understood in terms of the phenomenon
of aliasing. Thus the z-transform of the sequence of samples of a time
waveform is representative of the Laplace transform of the original
waveform in a way which is well understood. The Laplace transform of a
train of impulses repeats its values taken in a horizontal strip of the
s-plane of width 2ir/T in every other strip parallel to it. The z-transform
maps each such strip into the entire z-plane or, conversely, the entire
z-plane corresponds to any horizontal strip of the s-plane, for example,
the region -co < a < a>, —r/T ^ a> < ir/T, where s = a -+- jo.
In the same correspondence, the ;'co axis of the s-plane, along which
we generally equate the Laplace transform with the Fourier transform,
is the unit circle in the z-plane; the origin of the s-plane corresponds to
z ■* 1. The interior of the z-plane unit circle corresponds to the left
CHIRP Z-TRANSFORM
1251
half of the s-plane; the exterior corresponds to the right half plane.
Straight lines in the s-plane correspond to circles or spirals in the z-plane.
Figure 1 shows the correspondence of a contour in the s-plane to a
contour in the z-plane. To evaluate the Laplace transform of the im-
pulse train along the linear contour is to evaluate the z-transform of
the sequence along the spiral contour.
Values of the z-transform are usually computed along the path
corresponding to the ;« axis, namely the unit circle. This gives the
discrete equivalent of the Fourier transform and has many applications
including the estimation of spectra, filtering, interpolation, and corre-
lation. The applications of computing z-transforms off the unit circle
are fewer, but one is presented in this paper, namely the enhancement of
spectral resonances in systems for which one has some foreknowledge
of the locations of poles and zeros.
Just as we can only compute equation (2) for a finite set of samples,
so we can only compute equation (2) at a finite number of points,
say z fc :
x k - x(z k ) = £ «
f k \
(4)
The special case which has received the most attention is the set of
points equally spaced around the unit circle,
->■■ - ( ' X I' U ' j^*
k = 0, 1
,N - I
(5)
z-PLANE
S-PLANE
Fig. 1 — The correspondence of a z-plane contour to an s-plane contour through
the relation z = e' T .
1252 THE BELL SYSTEM TECHNICAL JOURNAL, MAY-JUNE 1969
for which
X k = E^expl-jf^nft] , k = 0, 1, ••• ,N - 1. (6)
n-(l
Equation (6) is called the discrete Fourier transform. The reader may
easily verify that, in equation (5), other values of k merely repeat
the same N values of z k , which are the iVth roots of unity. The discrete
Fourier transform has assumed considerable importance, partly be-
cause of its nice properties but mainly because since 1965 it has be-
come widely known that the computation of equation (6) can be
achieved, not in the N 2 complex multiplications and additions called
for by direct application of equation (6), but in something of the
order of N log2 iV operations if N is a power of two, or N 2 m* opera-
tions if the integers ?w< are the prime factors of N. Any algorithm
which accomplishes this is called a fast Fourier transform. Much of
the importance of the fast Fourier transform is that the discrete
Fourier transform may be used as a stepping stone to computing
lagged products such as convolutions, autocorrelations, and cross cor-
relations more rapidly than before. 3 - 4 The discrete Fourier transform
has some limitations which can be eliminated using the chirp 2-trans-
form algorithm which we describe. We also investigate the computa-
tion of the 2-transform on a more general contour, of the form
z k = AW~ k , k = 0, 1, • • • , M - 1 (7a)
where M is an arbitrary integer and both A and W are arbitrary
complex numbers of the form
A = A exp O'27r0 o ) (7b)
and
W = "Pr o exp0"27r^). (7c)
(See Fig. 2.) The case A = 1, M = N, and W = exp (—j2ir/N) corre-
sponds to the discrete Fourier transform. The general z-plane contour
begins with the point z = A and, depending on the value of W, spirals
in or out with respect to the origin. If W = 1, the contour is an arc of a
circle. The angular spacing of the samples is 2ir(p . The equivalent
s-plane contour begins with the point
s = o + jw = ^ In A (8)
CHIRP Z -TRANSFORM
1253
S-PLANE
Acr=y,lnW
CT = y/nA
Fig. 2 — An illustration of the independent parameters of the chirp z-transform
algorithm. The upper figure shows how the z-transform is evaluated on a spiral
contour starting at the point z = A. The lower figure shows the corresponding
straight line contour and independent parameters in the s-plane.
and the general point on the s-plane contour is
s k = s + 7c(Acr + ; Aco) = ^(ln A - k In W),
k = 0, 1, ••• ,M - 1. (9)
Since A and W are arbitrary complex numbers we see that the points
s fc lie on an arbitrary straight line segment of arbitrary length and
sampling density. Clearly the contour indicated in equation (7a) is
not the most general contour, but it is considerably more general
than that for which the discrete Fourier transform applies. In Fig. 2,
an example of this more general contour is shown in both the z-plane
and the s-plane.
To compute the z-transform along this more general contour would
seem to require NM multiplications and additions since the special
symmetries of exp (j2irk/N) which are exploited in the derivation of
the fast Fourier transform are absent in the more general case. How-
ever, we shall see that by using the sequence TF nV2 in various roles
we can apply the fast Fourier transform to the computation of the
z-transform along the contour of equation (7a). Since for W = 1, the
sequence W n ' /2 is a complex sinusoid of linearly increasing frequency,
and since a similar waveform used in some radar systems has the
picturesque name "chirp", we call the algorithm we are about to present
the chirp z-transform. Since this transform permits computing the
z-transform on a more general contour than the fast Fourier transform
1254 THE BELL SYSTEM TECHNICAL JOURNAL, MAY-JUNE 1969
permits, it is more flexible than the fast Fourier transform, although it
is also considerably slower. The additional freedoms offered by the
chirp 2- transform include:
(i) the number of time samples does not have to equal the number
of samples of the z-transform.
(it) Neither M nor N need be a composite integer.
(Hi) The angular spacing of the z k is arbitrary.
(iv) The contour need not be a circle but can spiral in or out with
respect to the origin. In addition, the point z„ is arbitrary, but this is
also the case with the fast Fourier transform if the samples x n are multi-
plied by z~ n before transforming.
II. DERIVATION OF THE CHIRP Z-TRANSFORM
Along the contour of equation (7a), equation (4) becomes
X h - S x n A' n W nk , h - 0, 1, ■ •• , M - 1 (10)
n-0
which, at first appearance, seems to require NM complex multiplica-
tions and additions, as we have already observed. But, let us use
Bluestein's ingenious substitution 5
nk = »' + *'-(*- »)' (11)
for the exponent of W in equation (10). This produces an apparently
more unwieldly expression
x k = T l x n A- n w (n,/2) w ik ' /2) w-«- n) ' /2 ,
n-0
k = 0, 1, ••• ,M - 1 (12)
but in fact equation (12) can be thought of as a three step process
consisting of: (i) forming a new sequence y n by weighting the x n ac-
cording to the equation
y n = x n A-'W" /a , n - 0, 1, • • • , N - I, (13)
(u) convolving y„ with the sequence v n defined as
Vn = w~ n ' /2 (14)
to give a sequence g k
g»
= Ep.- fc = 0,l, ■•• ,M-1, (15)
CHIRP Z -TRANSFORM
1255
and (Hi) multiplying g k by W*' /a to give X k
X k = g k W
k = 0, 1, • ■ • , M - 1.
(16)
The three step process is illustrated in Fig. 3. Steps i and Hi require
N and M multiplications respectively; step ii is convolution which may
be computed by the high speed technique disclosed by Stockham, based
on the use of the fast Fourier transform. 3 Step ii is the major part of
the computational effort and requires a time roughly proportional to
(N + M) log (N + M).
Bluestein used the substitution of equation (11) to convert a dis-
crete Fourier transform to a convolution as in Fig. 3. The linear sys-
tem to which the convolution is equivalent can be called a chirp filter
which is, in fact, also sometimes used to resolve a spectrum. Blue-
stein showed that for N a perfect square, the chirp filter could be
synthesized recursively with N^ multipliers and the computation of a
discrete Fourier transform could then be proportional to N 3/2 (see
Ref. 5) .
The flexibility and speed of the chirp 2-transform algorithm are
related to the flexibility and speed of the method of high-speed con-
volution using the fast Fourier transform. Recall that the product of
the discrete Fourier transforms of two sequences is the discrete Fourier
transform of the circular convolution of the two sequences; therefore,
a circular convolution is computable as two discrete Fourier trans-
forms, the multiplication of two arrays of complex numbers, and an
inverse discrete Fourier transform, which can also be computed by
the fast Fourier transform. Ordinary convolutions can be computed
as circular convolutions by appending zeroes to the end of one or
both sequences so that the correct numerical answers for the ordinary
convolution can result from a circular convolution.
We now summarize the details of the chirp z-transform algorithm
on the assumption that an already existing fast Fourier transform
n = o,i,...,N-i
<*>
yn
,nZ/2
HIGH SPEED
CONVOLUTION
,-r>2/2
y n ®w~
9k
x k
k = o,i,2,...M-i
W
k 2 /z
Fig. 3 — An illustration of the steps involved in computing values of the
z-transform using the chirp z-transform algorithm.
1256 THE BELL SYSTEM TECHNICAL JOURNAL, MAY-JUNE 1969
program (or special purpose machine) is available to compute discrete
and inverse discrete Fourier transforms.
We begin with a waveform in the form of N samples x n and we
seek M samples of X k where A and W have also been chosen:
(i) We choose L, the smallest integer greater than or equal to N +
M — 1 which is also compatible with our high speed fast Fourier
transform program. For most users this will mean L is a power of two.
Notice that while many fast Fourier transform programs will work
for arbitrary L, they are not equally efficient for all L. At the very
least, L should be highly composite.
(ii) We form an L point sequence y n from x n by the equation
\A- n W n ' /2 x n n = 0, 1, 2, ■ • • , N - 1 n ~
y n = \ • (17)
lo n = N,N + 1, ••• ,L - 1
(Hi) We compute the L point discrete Fourier transform of y n by
the fast Fourier transform, calling it Y r , r = 0, 1, . . . , L — 1.
(iv) We define an L point sequence v n by the relation
far*'" g n g M - 1
Vn = \w-' L - n) ' n L-N + l^n<L. (18)
arb itrary other n , if any
If L is exactly equal to M + N — 1, the region in which v n is arbitrary
will not exist. If the region does exist an obvious possibility is to in-
crease M, the desired number of points of the z-transform we com-
pute, until the region does not exist.
Notice that v n could be cut into two with a cut between n = M — 1
and n = L— N-\-V, if the two pieces were abutted differently, the result-
ing sequence would be a slice out of the indefinite length sequence
W~ n ' /2 . This is illustrated in Fig. 4. The sequence v n is defined the way
it is in order to force the circular convolution to give us the desired
numerical results of an ordinary convolution.
(v) We compute the discrete Fourier transform of v„ and call it V r ,
r = Q,l,...,L-l.
(vi) We multiply V r and Y r point by point, giving G> :
G T = V T Y r , r = 0, 1, • • • , L - 1.
(vii) We compute the L point inverse discrete Fourier transform
0fc,of G>.
Xn
(a)
(b)
(c)
yn
N-l
Y r
N-l
*W
-n2/2
L-i r
(d)
r-e- ARBITRARY-
M-!
L-N + i L-t n
♦ Vr
(f)
o:
(h)
(L)
L-l
♦ 9k
L-i r
NOT USED
M-1
L-I k
fX k
M-1
Fig. 4 — Schematic representation of the various sequences involved in the
chirp z-transform algorithm: (a) input sequence x„ with N values, (b) weighted
input sequence y n = A _ "TT'" 2 n x n ■ (c) discrete Fourier transform of y„. (d) values of
the indefinite sequence W~ ni n . (e) sequence v n formed appropriately from segments
of W~" 2 n . (f) discrete Fourier transform of v n . (g) product G> = Y r -V r . (h) inverse
discrete" Fourier transform of G> . (i) desired M values of the z-transform.
1258 THE BELL SYSTEM TECHNICAL JOURNAL, MAY-JUNE 1969
(viii) We multiply g k by W k%/2 to give us the desired X k :
X k = W k ' /2 g k , fc = 0, 1,2, ... ,M- 1.
The ^ for fc ^ M are discarded.
Figure 4 represents typical waveforms (magnitudes shown, phase
omitted) involved in each step of the process.
III. FINE POINTS OF THE COMPUTATION
3.1 Operation Count and Timing Considerations
An operation count can be made, roughly, from the eight steps just
presented :
(i) We assume that step i, that is, choosing L, is a negligible opera-
tion.
(ii) Forming y n from x n requires N complex multiplications, not
counting the generation of the constants A~ n W n%n . The constants may
be prestored, computed as needed, or generated recursively as needed.
The recursive computation would require two complex multiplications
per point.
(Hi) An L point discrete Fourier transform requires a time /c F ft^
\0g2L for L a power of two, and a very simple fast Fourier transform
program. More complicated (but faster) programs have more com-
plicated computing time formulas.
(iv), (v) The value of v n is computed for either M or N points, which-
ever is greater. The symmetry in W~ n ' /2 permits the other values of v n
to be obtained without computation. Again, v n can be computed re-
cursively. The fast Fourier transform takes the same time as that in
step Hi. If the same contour is used for many sets of data, V r need only
be computed once, and stored.
(vi) This step requires L complex multiplications.
(vii) This is another fast Fourier transform and requires the same
time as step Hi.
(viii) This step requires M complex multiplications.
As the number of samples of x n or X k grows large, the computation
time for the chirp ^-transform grows asymptotically as something
proportional to L log 2 L. This is the same sort of asymptotic depend-
ence of the fast Fourier transform, but the constant of proportionality
is bigger for the chirp 3-transform because two or three fast Fourier
transforms are required instead of one, and because L is greater than
N or M. Still, the chirp z-transform is faster than the direct com-
CHIRP Z -TRANSFORM 1259
putation of equation (10) even for relatively modest values of M
and N of the order of 50.
3.2 Reduction in Storage
The chirp z-transform can be put into a more useful form for com-
putation by redefining the substitution of equation (11) to read
. (n - N ) 2 + k 2 - (fc - n + N.) 2 + 2N k
nk- 2
Equation (12) can now be rewritten as
X„ = W h ' /2 W"' k 2! Xn A- n W tn - I, ' ) ' /a W- tk " ,+lr ' i ' / \
n-0
The form of the new equation is similar to equation (12) in that the
input data x n are preweighted by a complex sequence (A" n TF (n ~ Ar ° >,/2 ),
convolved with a second sequence (W~ ln ~ N ' ) ' /2 ) i and postweighted by a
third sequence (W k ' /2 W Nok ) to compute the output sequence X k .
However, there are differences in the detailed procedures for realizing
the chirp z-transform. The input data x n can be thought of as having
been shifted by N samples to the left. For example, x„ is weighted by
W N °' /2 instead of W°. The region over which W~ n ' must be formed, in
order to obtain correct results from the convolution, is
-N + I + N ^ n ^ M - 1 + N .
By choosing N = (N — M)/2 it can be seen that the limits over which
W~ n ' /2 is evaluated are symmetric; that is, W~ n ' /2 is a symmetric func-
tion in both its real and imaginary parts. (Therefore, the transform of
W~" 1/2 is also symmetric in both its real and imaginary parts.) It can
be shown that by using this special value of N a only (L/2 + 1) points
of W~ n ' /2 need be calculated and stored, and these (L/2 -f- 1) complex
points can be transformed using an L/2 point transform*. Hence the
total storage required for the transform of W~ "" is L + 2 locations.
The only other modifications to the detailed procedures for evaluating
the chirp z-transform presented in Section II are:
(i) following the L point inverse discrete Fourier transform of step
vii, the data of array g k must be rotated to the left by N locations,
\ii) the weighting factor of the g k is W k ' /2 W Nck rather than W k ' /a .
* The technique for transforming two symmetric L point sequences using one
L/2 point fast fourier transform was demonstrated by J. Cooley at the fast
Fourier transform workshop, Arden House, Harriman, New York, October 1968.
The appendix summarizes this technique.
1260 THE BELL SYSTEM TECHNICAL JOURNAL, MAY-JUNE 1969
The additional factor W N ° k represents a data shift of N samples to
the right, thus compensating for the initial shift and keeping the
effective positions of the data invariant to the value of N used.
Now we can estimate the storage required to perform the chirp z-
transform. Assuming that the entire process is to take place in core,
storage is required for V r which takes L + 2 locations, for y n , which
takes 2L locations, and perhaps for some other quantities which we wish
to save, such as the input or values of W*' /a or A~*W*' /a .
3.3 Additional Considerations
Since the chirp z-transform permits M ?* N, it is possible that occa-
sions will arise where M ~S> N or N » M. In these cases, if the smaller
number is small enough, the direct method of equation (10) is called
for. However, if even the smaller number is large it may be appropriate
to use the methods of sectioning described by Stockham. 3 Either the
lap-save or lap-add methods may be used. Sectioning may also be used
when problems, too big to be handled in core memory, arise. We have not
actually encountered any of these problems and have not programmed
the chirp z-transform with provision for sectioning.
Since the contour for the chirp z-transform is a straight line seg-
ment in the s-plane, it is apparent that repeated application of the
chirp z-transform can compute the z-transform along a contour which
is piecewise spiral in the z-plane or piecewise linear in the s-plane.
Let us briefly consider the chirp z-transform algorithm for the case
of z k all on the unit circle. This means that the z-transform is like a
Fourier transform. Unlike the discrete Fourier transform, which by
definition gives N points of transform for N points of data, the chirp
z-transform does not require M = N. Furthermore the z k need not
stretch over the entire unit circle but can be equally spaced along an
arc. Let us assume, however, that we are really interested in comput-
ing the N point discrete Fourier transform of N data points. Still the
chirp z-transform permits us to choose any value of N, highly com-
posite, somewhat composite, or even prime, without strongly affecting
the computation time. An important application of the chirp z-trans-
form may be computing discrete Fourier transforms when N is not a
power of two and when the program or special purpose device avail-
able for computing discrete Fourier transforms by fast Fourier trans-
form is limited to when N is a power of two.
There is also no reason why the chirp z-transform cannot be ex-
tended to the case of transforms in two or more dimensions with simi-
lar considerations. The two dimensional discrete Fourier transform
CHIRP Z-TRANSFORM 1261
becomes a two dimensional convolution which can be computed by-
fast Fourier transform techniques.
Caution: For the ordinary fast Fourier transform the starting point
of the contour is still arbitrary; merely multiply the waveform x n by
A~ n before using the fast Fourier transform and the first point on the
contour is effectively moved from z = 1 to z = A. However, the con-
tour is still restricted to a circle concentric with the origin. The angular
spacing of z k for the fast Fourier transform can also be controlled to
some extent by appending zeros to the end of x n before computing the
discrete Fourier transform (to decrease the angular spacing of the
z k ) or by choosing only P of the N points x n and adding together all
the x n for which the n are congruent modulo P; that is, wrapping the
waveform around a cylinder and adding together the pieces which
overlap (to increase the angular spacing) .
IV. APPLICATIONS OF THE ALGORITHM
Because of its flexibility, the chirp 2-transform algorithm discussed
in the Section III has many potential applications.
4.1 Enhancement of Poles
One advantage of the chirp 2-transform algorithm over the fast
Fourier transform is its ability to evaluate the 2-transform at points
both inside and outside the unit circle. This is important in the investi-
gation of systems whose transfer functions can be represented as ratios
of polynomials in z; that is, in finding poles and zeros of a linear sys-
tem. By evaluating the transform off the unit circle, one can make the
contour pass closer to the poles and zeros of the system, thus effectively
reducing the bandwidths and sharpening the transfer function.
For example, a five-pole system was simulated at a 10 kHz sam-
pling frequency. The poles were located at center frequencies of 270,
2290, 3010, 3500 and 4500 Hz with bandwidths of 30, 50, 60, 87 and
140 Hz, respectively. The z-plane pole positions are shown in Fig. 5.
(Those familiar with speech will recognize these numbers as resonance
positions appropriate for the vowel i in the word beet.) The input to
the system was a periodic impulse train of period 100 samples; that is,
100 pulses per second. Impulse invariant techniques were used to
simulate the system. 6 The 2-transform of one period of steady state
data (100 samples) was evaluated on two spirals outside the unit
circle, one on the unit circle, and two spirals inside the unit circle.
Figure 6 shows the five contours as they would appear in the s-plane
and the s-plane pole positions. The contours are seen to be equi-
1262 THE BELL SYSTEM TECHNICAL JOURNAL, MAY-JUNE 1969
Z- PLANE
Fig. 5 — Representation of the z-plane locations of the poles of the linear sys-
tem simulated in the text.
angularly spaced. The five sets of magnitude curves are shown in Fig.
7. The transform was evaluated at 50 equally spaced points from
to 4900 Hz, corresponding to <p = — 1/100. The sharpening of the
magnitude response in the region of the poles is quite pronounced.
Figure 6 indicates that contour 5 should be near optimum since it
intersects three of the poles.
This example is a somewhat idealized case in that spectral samples
were taken every 100 Hz; that is, at the harmonics of the funda-
mental frequency. Figure 8 shows the case for spectral data taken
every 25 Hz on contour 5 of Fig. 6, along with the case where the
spectral resolution is the same as shown in Fig. 7. This figure places
in evidence the true nature of the z-transform of a finite number of
samples. It is clear from equation (2) that X(z) has no poles any-
where in the z-plane except at z = 0. There are instead N—l zeros
which manifest themselves in the ripples seen in the upper curve of
Fig. 8. In many cases the poles of the original system which generated
the samples are still in evidence because the zeros tend to be arrayed
at approximately equal angular increments except at the locations of
the original poles. Hence a pole usually manifests itself by an absence
of zeros in the vicinity of that pole in the z-plane. Zeros of transmis-
sion are often masked by these effects when only a finite number of
samples are transformed. Examples of this effect are given after equa-
tion (23).
It is of interest to examine the ability of the chirp z-transform
algorithm to determine the bandwidth of a pole as well as its center
CHIRP Z -TRANSFORM
1263
frequency. To investigate this point, synthetic samples were generated
with the bandwidth of the lowest pole (Bx) variable from 10 to 320
Hz by factors of 2; all other bandwidths and center frequencies were
held at values used in the previous example. Again a fixed 100 pulse
per second source excited each of the systems. Figure 9 illustrates the
six sets of poles and three contours used in the investigation. Contour
3 extends into the right half-plane (spiral outside the unit circle) and
is only close to the lowest pole. Contour 2 corresponds to the unit circle
in the z-plane (that is, the discrete Fourier transform). Contour 1 is
an appropriate left-half plane contour (spiral inside the unit circle)
used for investigating this system. The resulting set of 18 magnitude
curves (six different sets of poles and three contours) are shown in
Fig. 10. The rows of Fig. 10 show magnitude curves with a fixed band-
width and variable contour, whereas the columns show curves on the
same contour with variable bandwidths. There are 801 points plotted
in each curve in the range to 5,000 Hz.
Looking down any column it is seen that as B± increases, the level
of the first resonance decreases steadily. The variation in fine spectral
detail resulting from the distribution of zeros in the neighborhood of
the poles of the original system is seen clearly in column 1. For ex-
4000
X
z 3000 \-
O 2000
s
1000
LEFT HALF-PLANE
-150
-50 50 100
FREQUENCY IN HZ
Fig. 6 — The s-plane locations of the poles of Fig. 5 and five contours for
evaluation of the z-transform.
1264
THE BELL SYSTEM TECHNICAL JOURNAL, MAY-JUNE 1909
1000 2000 3000 4000 5000
FREQUENCY IN HZ
Fig. 7 — Magnitude curves corresponding to evaluation of the z-transform on
the five contours of Fig. 6.
1000 2000 3000 4000 5000
FREQUENCY IN HZ
Fig. 8 — A comparison of high resolution and low resolution evaluations of
the z-transform. The spacing of points is 25 Hz in the upper curve and 100 Hz
in the lower curve.
CHIRP Z-TRANSFORM
1265
4000
a 2000 -
-320
-50
FREQUENCY IN HZ
150
Fig. 9 — The pole locations and contours used to investigate the possibility of
bandwidth determination using the chirp z-transform.
ample, the fifth resonance at 4,500 Hz is difficult to find in the upper
plots and almost missing in the lower plots, because of the presence of
a zero at the pole position. Furthermore, the frequency at which the
magnitude is minimum, that is, the closest zero to contour 1 shifts
from 2,500 to 2,700 to 800 to 1,100 Hz as bandwidth increases.
The plots in columns 2 and 3 show little or no variation from about
2,000 to 5,000 Hz where the appropriate contours are generally far
away from the zeros of the distributions. The resonance at 4,500 Hz
is always easy to locate on these plots, thus indicating the desirability
of both detailed close-up examination of the transform (as on con-
tour 1) and less detailed, further away looks at the magnitude curve
(as on contours 2 and 3). The magnitude curves in the regions to
2,000 Hz are still fairly sensitive to the exact zero distribution for
contour 2, and slightly sensitive for contour 3. It would appear from
Fig. 10 that there are cases when bandwidth can be determined either
from the width or the magnitude of the resonance. Further investiga-
tion is necessary before quantitative techniques for determining band-
widths are available.
The choice of the optimum contours is highly dependent on the
locations of the poles of the original system. In general there is no
1266
THE BELL SYSTEM TECHNICAL JOURNAL, MAY-JUNE 1969
I
N
X
l\l
II
II
DO
CO
snaanaa ni aaruiNsvw
CHIRP Z -TRANSFORM
1267
single contour on which all the poles are located since these contours
are essentially lines of constant Q = (center frequency) /(bandwidth).
Hence the choice of contour is highly dependent on which of the sys-
tem poles is of most interest in the particular problem. However, some
interesting observations can be made from studying magnitude curves
for systems whose poles are constant Q poles. Such a system was
simulated by keeping the pole center frequencies at the values used
previously and setting the Q of each of the poles to 20. A 100 pulse
per second impulse train was again used to excite the system and one
period of steady state data was analyzed along contours 1 to 7 shown
in Fig. 11. The pole positions of the system are shown in this figure
and are seen to coincide with contour 6 exactly. The magnitude curves
for these seven contours are shown in Fig. 12. These are high resolution
spectra containing 801 points from to 5,000 Hz. Notice that both
magnitude curves 5 and 6 accentuate the poles equally except in the
region of the fifth pole where curve 5 appears slightly better than
curve 6. The fact that this occurs is not surprising in view of the fact
that the pole is really manifested by an absence of a zero in an array
of zeros of approximately the same magnitude and angular spacing.
Another anomaly which can be attributed to the way in which the
4000 -
2 2000 -
-200 -150 -100
FREQUENCY IN HZ
Fig. 11 — The s-plane locations of constant Q poles and contours on which the
z-transform was evaluated.
1268
THE BELL SYSTEM TECHNICAL JOURNAL, MAY-JUNE 1969
90 -,
70
50
30 h
105
85
65
45
1000 2000 3000 4000
FREQUENCY IN HZ
Fig. 12 — Magnitude curves for contours of Fig. 11 (constant Q poles).
zeros are distributed is evidenced by comparing curves 5 and 7. Based
on the relative positions of the two contours with respect to the poles
we would expect the magnitude curves for these contours to be identi-
cal, but the comparison shows that this is not the case in actual com-
putation. This is the result of the fact that the zero distribution is not
exactly symmetric so that contours which pass very close to the zeros
may look considerably different from one another.
A final point of interest in Fig. 12 is the linear component in the
last three curves which dominates at high frequencies. This effect is
also shown in Fig. 14. Figure 13 shows the five contours used in ob-
CHIRP Z -TRANSFORM
1209
taining the log magnitude plots in Fig. 14. It is clear that when the
contour passes inside the original pole locations (and therefore inside
the array of zeros in the 2-plane) , the log magnitude function exhibits
a definite linear component. This effect is easily explained when X(z)
is written in the form
X(z) - Dz-' N ~ l) 11(1 -a7 l z).
(19)
where the a r 's are the zeros of X(z). If we evaluate equation (19) at
z = W~ k exp (—j2TT(p k) we obtain for the magnitude
\X k \= \D\ WV' l} II 111 - a; l W: h exp i-frrJk)]\- (20)
For plotting in dB we define
20 log,„ | X k | = 20 log.o \D\ + 20(N - l)k log 10 W
+ Z 20 log 10 |[1 - a: x W: k exp (-j2ir^.fc)]|. (21)
r-l
In the examples we have shown, almost all of the zeros a r have mag-
nitudes slightly less than 1. Thus for contours inside these zeros (cor-
9, 2000
-600 -500 -400 -300 -200
FREQUENCY IN HZ
Fig. 13 — Contours and pole locations used to study the effect of passing in-
side the pole locations.
1270
THE BELL SYSTEM TECHNICAL JOURNAL, MAY-JUNE 1969
85
A A A .
70
55
A MV/V\
i i i i
110
90
70
^ — r^ i i i
140
no
80
60
■ / ^— ■ -r i i i
240
^^^
180
^ — ^^
120
^^~^
60
— -<^ i i i
280
200
120
40
•* i i i i
NO. 2
NO. 3
NO. 4
1000 2000 3000 4000 5000
FREQUENCY IN HZ
Fig. 14 — Magnitude curves for the contours of Fig. 13 showing a large linear
component resulting from the N — 1 poles of X(z) at z = 0.
responding to contours 2 through 5 in Fig. 13), W is greater than 1.
Thus each term in the sum on the right side of equation (21) tends to
decrease as k gets larger. In contrast, the second term on the right
side of equation (21) represents a linear component with slope equal
to20 (iV— 1) log 10 F .
In Fig. 14, the rath curve corresponds to a value of W = e 47rm/10 '° 00 ,
m = 1, 2, 3, 4, 5. The value of <p is - 1/100, N is 100, and the sam-
pling rate is 10 kHz. Thus the frequency going from to 5 kHz cor-
responds to k going from to 50. For example, in the fifth curve W =
e 2 0w /io.ooo anc j t h e s i p e should be
20(tf - 1) log.. W. = ^y" - 5.4. (22)
Thus the total dB change in going from to 5 kHz should be on the
order of 50(5.4) = 270 dB. In Fig. 15a we show this case again. In
Fig. 15b we show the result of evaluating
2; w.r = 2 wr " I, Ki
(23)
»— (AT-l)
CHIRP Z-TRANSFORM
1271
using the same value of W and <p . Notice that this should remove the
second term in equation (21) leaving the other terms unaffected. This
observation is substantiated by Fig. 15b since the value at 5 kHz is
very nearly 270 dB less than in Fig. 15a. Notice also that some of the
resonances are still in evidence although not as clearly defined because
the contour is passed relatively far inside the zeros of X(z).
Another interesting question was investigated using this technique.
As shown above, when the response of a linear system is truncated by
repetitively pulsing the system and transforming a finite number of
samples, the z-transform has only zeros except for poles at z = 0.
However, the poles of the original system function can still be located
in magnitude response curves by the absence of zeros in the ap-
propriate regions. The question arises about what happens when the
system function contains zeros. Suppose h (nT) is the impulse response
of a linear system with both poles and zeros in its z-transform H(z).
Since H(z) has poles, h(nT) will be an infinite sequence. There is
clearly no reason to expect that the transform of only N of these
samples will have zeros at the same location as the zeros of H(z).
However, the system zeros can be expected to have an effect on the
distribution of zeros of the truncated z-transform.
To illustrate this point, the system response used in the above ex-
amples was modified by passing the output waveform through a sys-
tem whose transfer function consisted of a complex conjugate zero
pair. A periodic 100 pulse per second source was again used to excite
the system and one period of steady state data was analyzed. The
system pole-zero pattern and the contours of analysis are shown in
I0OO 2000 3000 4000 5000
FREQUENCY IN HZ
Fig. 15 — Magnitude curves obtained by evaluation of the z-transform on
contour 5 of Fig. 13: (a) with the effect of the N — 1 poles at z = 0; (b) with
the N — 1 poles removed by shifting the sequence x» by N positions to the
left.
1272 THE BELL SYSTEM TECHNICAL JOURNAL, MAY-JUNE 1969
Fig. 16. In one simulation a zero was placed at point A (500 Hz, 12.5
Hz), and in a second simulation the zero was at point B, (2,500 Hz,
60 Hz). The analysis was made at 801 points from to 5,000 Hz along
these contours. The resulting magnitude curves, along with a set of low
resolution curves where the magnitude was computed every 100 Hz
from to 5,000 Hz, are shown in Figs. 17 and 18. The data of Fig. 17
are for the case where the transmission zero was at 500 Hz whereas
for Fig. 18 the zero was at 2,500 Hz.
The high resolution data of Fig. 17 show no strong indication of
the transmission zero; whereas the transmission poles are still very
much in evidence. The low resolution data (evaluated at harmonics
of the source) does indicate the presence of a zero along contour 1,
but along the other contours the case is not so clear. The most unusual
observation is that along contour 3, the contour closest to both the
transmission zero and the poles, there is little or no indication of the
zero; whereas the poles are still strongly in evidence. Along contour
4, at the high frequency end, there is noise in the magnitude spectrum.
The source of this noise is discussed in Section V.
The indications from Fig. 17 are that a transmission zero can be
more easily located on contours which are far from the zero than on
5000
4000 -
9, 2000 -
1000 -
-120 -90 -60
FREQUENCY IN HZ
Fig. 16 — The s-plane locations of poles and zeros (at A and B) and contours
used in studying the effect of zeroes on the magnitude curves.
CHIRP Z-TRANSFORM
1273
HIGH RESOLUTION DATA
(801 POINTS)
190
160
130
100
" 1 1 1 1
190
A 7 Nw/ V_w^w
160
130
^
fvvttftnr
i • i i i
190
160
w-v-
130
*y
««*^
1 1 1 1
280
240
200
160
IPO
-~"l 1 1 !
LOW RESOLUTION DATA
(51 POINTS)
1000 2000 3000 4000 5000
FREQUENCY IN HZ
NO.?
N0.4
5000
Fig. 17 — Magnitude curves for a zero at 500 Hz (position A in Fig. 16).
contours which traverse it. Furthermore it is much easier to locate on
a low resolution spectrum than on a high resolution spectrum. Hence
zeros, unlike poles, are not generally easy to locate from spectra.
The zero of Fig. 17 was at 500 Hz and in a region where the high
resolution spectra displayed a large amount of ripple from the trunca-
tion zeros of the data. Figure 18 shows similar magnitude curves for
the zero at 2,500 Hz, a region with much less ripple in the spectrum.
The magnitude response curves show effects entirely similar to those
of Fig. 17. The zero is most easily beatable for contour 1, the stand-
ard fast Fourier transform. In contour 3, which again passes through
the zero, there is no indication of the zero. Also the low resolution data
tends to show the zero better than the high resolution data. One im-
portant implication of these results is that one could not use these
techniques to accurately find the position of complex transmission
zeros. In many cases it would be difficult to differentiate between dips
in the spectrum between poles and dips caused by complex zeros, thus
indicating the difficulty of locating even the center frequency of a
zero.
The chirp z-transform algorithm has been applied to the spectral
analysis of speech in order to aid in automatic detection of the time
varying resonances (poles or formants) of speech. Voiced speech can
1274
THE BELL SYSTEM TECHNICAL JOURNAL, MAY-JUNE 1969
HIGH RESOLUTION DATA
(801 POINTS)
LOW RESOLUTION DATA
(51 POINTS)
N0.1
N0.2
NO.3
N0.4
1000 2000 3000 4000 5000 1000
FREQUENCY IN HZ
2000 3000 4000 5000
Fig. 18 — Magnitude curves for a zero at 2,500 Hz (position B in Fig. 16).
be modelled as the convolution of a source waveform with a vocal
tract impulse response. The vocal tract impulse response is essentially
a sum of damped exponentials, each exponential corresponding to a
mode or pole of the vocal tract transfer function. It is of interest to
speech researchers to detect these time varying resonances. The chirp
2-transform algorithm has been applied to individual periods of voiced
speech with a high degree of success. Figure 19 shows the result of ap-
plying the chirp 2-transform algorithm along the two contours shown
at the upper left of the figure, to a period of voiced speech. The upper
contour corresponds to the standard fast Fourier transform contour;
the lower to a suitably chosen spiral contour. The magnitude function
along the upper contour indicates a single wide peak in the region
2,000 to 2,500 Hz, whereas the magnitude along the lower contour
shows two isolated peaks in this region corresponding to the physical
knowledge that there actually are supposed to be two peaks in this
region. Variations on the chirp 2-transform algorithm for spectral
analyses of speech have been studied and will be reported on in a
subsequent paper. 7
4.2 High Resolution, Narrow Band Frequency Analysis
One very useful application of the chirp z-transform algorithm is the
ability to efficiently evaluate high resolution, narrow frequency band
CHIRP Z-TRANSFORM
1275
spectra. Using standard fast Fourier transform techniques, in order to
achieve a frequency resolution of ^ AF, with a sampling frequency of
the data of l/T, requires N ^ 1/(T-AF) points. For very small AF,
this implies very large values of N. The crucial issue is that what is
often required is high resolution for a limited range of frequencies and
low resolution for the remainder of the spectrum. An example of such a
circumstance is the design of band-pass or low-pass filters. Usually
what is desired is a microscopic look at details of the frequency response
in the pass-band and only a gross look outside the pass-band.
The chirp 2-transform algorithm is extremely well suited for such
cases since it allows selection of initial frequency and frequency spac-
ing, independent of the number of time samples. Hence high resolu-
tion data over a narrow frequency range can be attained at low cost.
To illustrate these points, simple rectangular band-pass niters were
simulated by symmetrically truncating a delayed impulse response.
The impulse response used was
h(nT) = a sin [ir(F 2 - F,)(n -\- m)T]
•cos [t(F 2 + F^n - \ - m)T),
< n < 2ra
(24)
S-PLANE
3500 4000
Fig. 19 — Magnitude curves from evaluation of the z-transform of one period
of natural speech. The contour for the upper plot is the unit circle in the 2-plane
while the contour for the lower curve is a spiral inside the unit circle.
1276 THE BELL SYSTEM TECHNICAL JOURNAL, MAY-JUNE 1969
where
2ra = number of terms in the truncated impulse response
1/T = sampling frequency = 10,000 Hz
F x = lower cutoff frequency in Hz
F 2 = upper cutoff frequency in Hz.
Values for m of 100 and 500 were used with F x = 900 Hz and F 2 =
1,100 Hz. Figure 20 shows plots of equation (24) for these two cases. A
standard 1,600 point fast Fourier transform was calculated and the
magnitude response for m = 100 is shown in the upper half of Fig. 21.
In order to investigate the pass-band and transition region more care-
fully the chirp z-transform algorithm was used to give a 1.25 Hz
resolution over the band from 500 to 1,500 Hz. The contour used was
identical to the contour for the fast Fourier transform. The resulting
magnitude response curve is shown in the lower half of Fig. 21. To
achieve this high a resolution would have required an 8,000 point
fast Fourier transform, instead of the 1,000 point transforms actually
used. (Similar expansions of regions of the phase curve were made for
this filter but are not shown.)
Figure 22 shows similar effects for the case m = 500. The applica-
bility of the chirp z-transform algorithm for such frequency expan-
sions is a powerful tool for close examination of small frequency
bands, as well as for debugging implementations of digital filters.
For example, one could easily check if a desired filter met its design
specification of in-band ripple, transition ratio, and so on. 8
One situation where the chirp z-transform algorithm may be quite
useful is when we are confronted with an extremely long sequence for
which we desire a fine grained spectrum over a narrow band of fre-
quencies. Suppose we have a sequence of P samples and desire M
spectral samples where M <<C P. That is, we wish to evaluate
p-\
X k = £ x n A~ n W k , k = 0, 1, ••• ,M - 1. (25)
n-0
The sum in equation (25) can be broken up into r sums over N points
as follows
X k m 2 A-" V W
kqX
2 x n+QlV A- n W nk , k - 0, 1, ■ • • , M - 1
n-0 J
(26)
where rN ^ P. Each of the r sums in the brackets can be evaluated
using the chirp z-transform algorithm, requiring storage on the order of
CHIRP Z-TRANSFORM
1277
10000
4^4^*
300 400 500 600
SAMPLE NUMBER
700 800 900 1000
Fig. 20 — The impulse responses of simple band-pass filters.
3(N + M — 1) locations. In addition we require 2M locations in which
to accumulate the M complex values of the transform. Although 2
fast Fourier transforms and 2M complex multiplications are required
for each of the r transforms, it is quite possible that a saving in total
time may result from this method as opposed to evaluation of a P
point transform using auxiliary storage such as drum, disk or tape.
1278
THE BELL SYSTEM TECHNICAL JOURNAL, MAY-JUNE 1969
160
- A
100
120
Mgf^ ^mmrnmm-
80
\\\\W m]\\mm
pPTOfm
PfRfftam^
in
i 4 °
o
Q
i i i i i
i i
i
2000 3000
FREQUENCY IN HZ
900 1000 1100
FREQUENCY IN HZ
Fig. 21 — Frequency response curves for upper impulse response (200 samples)
in Fig. 20. Upper curve obtained with 1,600 point fast Fourier transform (resolu-
tion 6.25 Hz). Lower curve obtained with chirp z-transform algorithm (1.25 Hz
resolution) .
4.3 Time Interpolation or Sampling Rate Changing
The flexibility of the chirp ^-transform algorithm for obtaining high
resolution in frequency has been explained and illustrated in Section
4.2. A similar procedure applies to interpolation between samples of a
bandlimited time function using samples of the frequency spectrum. 9
In this section, we discuss how the discrete Fourier transform can be
used to perform interpolation on a set of samples and the advantages
and disadvantages of using the chirp 2-transform algorithm for this.
4.3.1 Bandlimited Interpolation Using the discrete Fourier transform
Assume that we have available N samples x(nT) , n = 0, 1, 2, . . . ,
N — 1, of a bandlimited waveform x{t). The sampling interval T is
assumed less than or equal to the Nyquist interval. The total time in-
terval spanned by these samples is therefore NT seconds. We wish to
obtain equally spaced samples of x(t) at a sampling interval T f , where
CHIRP Z -TRANSFORM
1279
T is less than or equal to the Nyquist interval. These samples are
denoted by x(mT') , m = 0, 1, . . . , N' - 1, where N'T' = NT. (Notice
that we are assuming N' is an integer. This assumption will be dropped
later.)
If all the samples x(nT) are available, the samples x{mT') can be
obtained from
x(mT') = £ x ( nT )
sin ^ (mT r - nT)
| (mT' - nT)
(27)
Thus the interpolation can be viewed as the result of convolving the
interpolation function [sin (7rt/T)]/[(ir/T) t] with the samples re (nT)
and then resampling with period T'. It is well known that convolution
may be done using the discrete Fourier transform, and we will show
how the resampling can also be affected by properly augmenting the
160
- n
m= soo
120
mmmdtW \\
80
Wnvfv "\
TWnpfyw^^
m
w 40
u
* n
i i
M ' II ' l ' | j 1 '
i i i i i i i
2000 3000
FREQUENCY IN HZ
5000
900 1000 1100
FREQUENCY IN HZ
Fig. 22 — Frequency response curves for lower impulse response (1,000 sam-
ples) of Fig. 20. Upper curve obtained with 1,600 point fast Fourier transform
(resolution 6.25 Hz). Lower curve obtained with chirp z-transform algorithm
(1.25 Hz resolution).
1280 THE BELL SYSTEM TECHNICAL JOURNAL, MAY-JUNE 1969
transform with zeros. Because the discrete Fourier transform uses only
a finite number of samples, we shall encounter errors similar (but not
identical) to using only N terms in equation (27) .
The discrete Fourier transform of the given samples is
XM = g x(nT) exp (- j || knTJ
k = 0, 1, ••• ,N - 1. (28)
(Notice that we have changed notation in this section in order to make
explicit the number of samples and the sampling period.) We define
X&k) = X N {k)H N (k), (29)
where H N (k) is the N point discrete Fourier transform of the inter-
polation function to be convolved with the samples x(nT). [Notice
that this convolution is equivalent to cyclic convolution of a periodic
impulse response h(nT) with the samples x{nT) .]
In order to change the sampling to period T', we split X/,{k) about
k = N/2 and expand (by inserting zeros) or contract (by discarding
zeros) the transform according to the following equations
X' N .(k) = xm
10 ^ k < N'/2
lo ^ k < N/2
N' <N
N' > N
(30a)
=
k = N'/2
N' < N
(30b)
= 2Aat(/c)
k = N/2
N' > N
(30c)
=
N/2 <k <N' - N/2
N' > N
(30d)
= \X' N (k -
- N' + N) k = N' -
N/2
N' < N
(30e)
= X' N (k -
N , + J N>/2 <k<N>
IN' - N/2 < k < N'
N' <N^
N' > N
(30f)
Equations (30b) , (30c) , and (30e) are required only when N' and N
are even integers and equation (30d) is required only when N' > N.
The N' point inverse discrete Fourier transform of X' N , (k) is defined
to be
x'{mT) = ~ g XUk) exp (; ^ kmT') ,
m = 0, 1, ••• ,N f - 1. (31)
CHIRP Z-TRANSFORM 1281
For example, if N is even and N' > N, we can show using equations
(28), (29), and (30) that
I AT/2 Tat-i / 9 ^
x'(mT') = p k S *»<*)[ § *(») ex P (-J ^f ^
■exp(;^7/fcmr), (32)
where m = 0, 1, • • • , N' — 1, and the terms corresponding to k = ±N/2
are understood to be multiplied by 1/2 since N is even. By interchanging
the order of summation and using the fact that N'T' = NT, we obtain
x'imT') = tt/ E x(nT)h(mT' - nT) (33)
*V „-0
where
h(mT' ~nT) = jj ^ H„(k) exp j f^ k(mT' - nT) J ■ (34)
Notice that equation (33) has the desired form for interpolation [see
equation (27)], however the values of x'imT') are clearly not exactly
equal to the desired interpolated values x{mT / ), This is so because
only N samples are used and because of the form of h (mT — nT) . As
an example, suppose
HM = fl -N/2 < k < 2V/2. (35)
li k= ±N/2
(This is equivalent to splitting X N {k) at k = N/2 and inserting N' —
N zeros between the two halves of the transform) . If we evaluate
equation (34) for this case, we obtain
sin ~ (mT' - nT)
h(mT' - nT) = (36)
Ntanj^(mT'-nT)
This function is plotted in Fig. 23a where = (mT' — nT)/T, and
N = 8. Clearly
sin ~ (mT' - nT)
h(mT' - nT) « (37)
| (m2" - nT)
1282 THE BELL SYSTEM TECHNICAL JOURNAL, MAY-JUNE 1969
(a)
N = 8
K
t ^^°~~-
(b)
1 I
1 1
1
1
T
2T
1 1
NT
I
t
1
2
N
n
N' m
Fig. 23 — An illustration of bandlimited interpolation using the discrete Fourier
transform: (a) the periodic function which is convolved with the original sam-
ples; (b) a bandlimited time function showing samples with spacing T; (c) How
the interpolated values xfimT') are formed.
when 7r(mT' — nT)/NT\s small, that is, in a region where h(mT' — nT)
is significantly different from zero. Figure 23b shows a segment of a
waveform x{t) and samples x(nT). In Fig. 23c, we have shown just
two of the terms in equation (33). This figure places in evidence the
nature of the interpolation which is performed. The errors are likely
to be greatest at either end of the segment, since the interpolated values
at one end depend on the samples at the other end in a way which is
not at all consistent with equation (27). The error caused by this
effect will be most significant in the regions ^ mT' < 2T and
CHIRP Z-TRANSFORM
1283
T(N — 2) ^ mT' ^ TN. The remainder of the values will have es-
sentially the same error associated with using only N terms in equa-
tion (27).
Notice that equation (36) is not the only interpolation function
which can be used. Other choices of H N (k) may lead to interpolation
functions which are in some sense more desirable. For example, Fig.
24 shows four different choices for H N (k) and their associated inter-
polation functions or impulse responses. (The impulse responses were
shifted modulo N' as an aid in plotting.) It can be seen from Fig. 24,
that removing the sharp cutoff in H N {k) greatly shortens the effective
duration of the impulse response, thus tending to minimize the end
effects discussed previously. Clearly the approximation to (sin wt)/vt
interpolation is not as good as equation (36) , but in many cases such
smoothing of the interpolated values may not be objectionable.
1.0
0.8
o.e
0.4
0.2
800
w 600
3 400
£ 200
-200
_ IMPULSE
RESPONSE
FREQUENCY
RESPONSE
1000 2000 3000 4000
FREQUENCY IN Hz
20 30 40
TIME IN SAMPLES
(a)
1.0
0.8
0.6
0.4
0.2
800
u 600
400
( 200
-200
FREQUENCY
RESPONSE
1000 2000 3000 4000 5000
FREQUENCY IN HZ
IMPULSE
RESPONSE
40 60 80 100 120 140 160
TIME IN SAMPLES
(b)
1.0
0.8
0.6
0.4
0.2
800
600
UJ
3 400
£ 200
-200
FREQUENCY
RESPONSE
1000 2000 3000 4000
FREQUENCY IN HZ
IMPULSE
RESPONSE
20 30 40
TIME IN SAMPLES
(O
i.o
0.8
0.6
0.4
0.2
800
600
J
3 400
C 200
-200
FREQUENCY
RESPONSE
1000 2000 3000 4000 5000
FREQUENCY IN HZ
IMPULSE
RESPONSE
20 40 60 80 100 120 140 160
TIME IN SAMPLES
(d)
Fig. 24 — A set of four simple frequency responses and corresponding im-
pulse responses which could be used for interpolation.
1284
THE BELL SYSTEM TECHNICAL JOURNAL, MAY-JUNE 1969
4.3.2 Computational Considerations in Bandlimited Interpolation
In Section 4.3.1 we discussed a method of bandlimited interpolation
based on the discrete Fourier transform. The operations involved are
summarized in Fig. 25. The sequence {X N {k) } may be evaluated using
the fast Fourier transform. In this case, N will be restricted to a value
compatible with the available fast Fourier transform routine, for ex-
ample, N would be highly composite. The transform is then multipled
by H N (k) and expanded or contracted according to equation (30).
Then we must compute the inverse discrete Fourier transform with N'
points. This can be done using the fast Fourier transform provided
that
(i) JV' = NT/1" is an integer and compatible with the available fast
Fourier transform routine.
(ii) Enough high speed storage (a minimum of N' locations) is
available.
[Notice that i applies for either N' > N or N' < N while ii will probably
not be a problem except when N' y> N.]
In many cases it may not be possible to meet one or both of the
above conditions; then the chirp z-transform algorithm can be very
useful. The N' point inverse discrete Fourier transform may be com-
puted using W = 1 and <p = + l/iV', where N' need not even be an
integer. Thus we can compute M interpolated values using
N/2
2
k—N/2
1 \LLz
x'(mT) = ^7 W m ' n E [X' N .(k)W" /2 ]W- (
-(m-*)V2
where X£.(k) is determined by equation (30) and
X' K .(-k) = X' N .(N' ~k), k = l,2,
,N/2.
(38)
(39)
x(nT)
DISCRETE
FOURIER
TRANSFORM
H N (k)
(Xhr
x N (k) v-^ X ' N (k)
n = o,i,...,N-i
N'
EXPAND
OR
CONTRACT
xW(k)
N'
INVERSE
DISCRETE
FOURIER
TRANSFORM
x'(mT')
m = o,i,...,N-i
Fig. 25 — Illustration of the steps involved in bandlimited interpolation using
the discrete Fourier transform.
CHIRP Z-TRANSFORM 1285
Assuming that the transform of W~ k ' /2 is available, equation (38)
can be evaluated using two L point fast Fourier transforms where L
is the smallest integer which is greater than N + M — 1 and which is
compatible with an available fast Fourier transform routine.
Alternatively we can evaluate
»W = W *■"■ 2! [Y'Ak)W" /2 ]W- l - k) ' /2 (40)
1 v ,. =
where
Y> N .(k) = \ X ' Ak) k = ° • (41)
[2Xit.(k) < k ^ N/2
It can be shown that
y'{mT') = x'{mT') - jx'(mT') (42)
where x'(mT') is the inverse discrete Fourier transform of
l' N .{k) = jsgn„.(k)-X' N .(k) (43)
and
k = 0, N'/2
sgn jV - (fc) - I 1 < k < N'/2. (44)
-1 N'/2 < k < N' - 1
From equation (41) and (44) it can be shown that x~'(mT') is an approxi-
mation to the Hilbert transform of x'{mT). In this case we require at
least {N/2 + M — 1) point transforms to compute M interpolated
values. This is at the expense of not being able to do two interpolations
at once as is possible with equation (38) (obtaining one interpolation
as a real output and one as an imaginary output); however we do obtain
an approximation to the Hilbert transform of x'iinT') which may be of
value in some applications.
If sufficient core storage is not available to compute an N' point
fast Fourier transform, we can compute the interpolated values in
sections and piece these sections together as is commonly done in high
speed convolution. The chirp ^-transform algorithm allows us to com-
pute as many as 2M interpolated points at a time, where M can be
chosen so that the fast Fourier transforms can be done using only core
storage. Probably the most significant advantage, though, is the
ability to efficiently interpolate to arbitrary sampling intervals.
128G
THE BELL SYSTEM TECHNICAL JOURNAL, MAY-JUNE 1969
As an example of the ideas discussed in this section, consider the
waveforms shown in Fig. 26. Figure 26a shows 500 samples of a speech
waveform where the sampling rate was 20 kHz. (T = 5 X 10~ B sec-
ond). The samples are connected by straight lines in the figure. Figure
26b shows the 500 samples of the waveform in (a) after filtering with
a nonrecursive filter of the type shown in Fig. 24, whose gain was zero
after 3.2 kHz. Figure 26c shows 160 samples of the result of a change
of sampling rate from 20 kHz to 6.4 kHz. The value of N was 700 and
N' = (6,400) (700)/20,000 = 224. It is difficult to judge quantitatively
from such a figure the accuracy of the interpolation. It does seem safe
to conclude that the error is not extreme. Our experience has been that
there is significant error only in the first and last few samples of the
N' output samples. Using the chirp z-transform algorithm, these "bad"
samples need not ever be computed, for example, only M "good" values
need be computed.
10 15
TIME IN MILLISECONDS
Fig. 26 — An example of interpolation for the purpose of changing the sampling
rate: (a) 500 samples of speech at 20 kHz sampling rate; (b) 500 samples of (a)
after low pass filtering to 3.2 kHz; (c) 160 samples of (a) after changing the
sampling rate to 6.4 kHz using the chirp z-transform. (In all cases the samples
are connected by straight lines).
CHIRP Z-TRANSFORM 1287
If one wishes to low-pass filter a waveform and then go to a lower
sampling rate, the filtering and interpolation can be combined if we
use a nonrecursive filter. That is, the discrete Fourier transform of the
filter impulse response can be simply combined with H N (k).
V. LIMITATIONS
Several times we have pointed out shortcomings of the chirp z-trans-
form algorithm. One limitation in using it to evaluate the z-transform
off the unit circle stems from the fact that we may be required to com-
pute Wt n ' /2 for large n. If W differs very much from 1.0, TP*" V2 can
become very large or very small when n becomes large. (We
require a large n when either M or N become large, since we
need to evaluate W*' /a for n in the range — N < n < M.) For example, if
W = e - 2 - 5/l0 - 000 « 0.999749, and n = 1,000, JP*" va = e il28 which
exceeds the single precision floating point capability of most computers
by a large amount. Hence the tails of the functions W ±n ' /2 can be
greatly in error, thus causing the tails of the convolution (the high fre-
quency terms) to be grossly inaccurate. The low frequency terms of the
convolution will also be slightly in error, but these errors generally are
negligible.
An example of this effect is shown in Fig. 27. The contour for the
five curves in this figure was held fixed (contour 5 in Fig. 6) and the
number of frequency points in the range to 5,000 Hz was increased
in steps of 2 from 50 to 800. Spectral samples are plotted every 100 Hz
for comparison. (This example was programmed using single-precision
floating-point arithmetic on a GE 635 computer with a 36 bit word
length.) It is seen that as the number of output points increases, errors
in the high frequency region become large and completely mask the
fifth resonance for the 800 point case. The effects of the inaccuracy in
W ±n ' /2 can also be seen at low frequencies. For example, the spectral
magnitude at Hz goes from about 120 dB to 134 dB as the number of
points goes from 50 to 800. These small errors generally do not affect the
gross spectral characteristics as seen in Fig. 27. The resonances are
easy to locate in all cases until the errors get exceedingly large. One
can push the maximum point limit higher than 800 (in this case) by
using double precision arithmetic.
The limitation of contour distance in or out from the unit circle is
again the result of computation of W ±n ' /2 . As W deviates significantly
from 1.0, the number of points for which W ±n ' /2 can be accurately
computed decreases. It is of importance to stress, however, that for
1288
THE BELL SYSTEM TECHNICAL JOURNAL, MAY-JUNE 1969
1000 2000 3000 4000
FREQUENCY IN HZ
Fig. 27 — A comparison of magnitude plots for varying number of points on the
same spiral contour. The fifth plot shows the effect of errors in evaluating W n2/J for
large n. (Points are plotted every 100 Hz in each curve to aid in comparison.)
W = 1 there is no limitation of this type since W ±n ' /2 is always of
magnitude 1.
The other main limitation of the chirp z-transform algorithm stems
from the fact that two L point fast Fourier transforms and one L/2
point fast Fourier transform must be evaluated where L is the smallest
convenient integer greater than N + M — 1 as previously mentioned.
We need one fast Fourier transform and 2L storage locations for the
transform of x n A~ n W n ' /2 ; one fast Fourier transform and L+2 storage
locations for the transform of W~ n ' /2 ; and one fast Fourier transform
for the inverse transform of the product of these two transforms. We
do not know a way of computing the transform of W~ n ' /a either re-
cursively or by a specific formula (except in some trivial cases.) Thus
we must compute this transform and store it in an extra L + 2 storage
locations. Of course, if many transforms are to be done with the same
value of L we need not compute the transform of W~ n ' /2 each time.
We can compute the quantities A~ n W n7/2 recursively, as they are
CHIRP Z -TRANSFORM 1289
needed, to save computation and storage. This is easily seen from the
fact that
A- b+ V» ,fl - (A-W^-WlflA- 1 . (45)
If we define
and
then
and
C n = A-W n ' /2 (46)
D n - WW+A- 1 (47)
D n+l = W-D n (48)
C n+1 = C n -D„. (49)
Setting A = 1 in equations (45) through (49) provides an algorithm
for the coefficients required for the output sequence. A similar recursion
formula can be obtained for generating the sequence A~ n W in
The user is cautioned that recursive computation of these coefficients
may be a major source of numerical error, especially when W ~ 1
or <p tt 0.
VI. SUMMARY
We give a computational algorithm for numerically evaluating the
2-transform of a sequence of A r time samples. This algorithm, the
chirp 2-transform algorithm, enables the evaluation of the z-trans-
form at M equiangularly spaced points on contours which spiral in or
out (circles being a special case) from an arbitrary starting point in
the z-plane. In the s-plane the equivalent contour is an arbitrary
straight line.
The chirp 2-transform algorithm has great flexibility in that neither
A r or M need be composite numbers; the output point spacing is
arbitrary ; the contour is fairly general and N need not be the same as
M. The flexibility of the chirp 2-transform algorithm comes from being
able to express the 2-transform on the above contours as a convolu-
tion, permitting the use of well-known high speed convolution tech-
niques to evaluate the convolution.
Applications of the chirp 2-transform algorithm include enhance-
ment of poles for use in spectral analysis, high resolution narrowband
1290 THE BELL SYSTEM TECHNICAL JOURNAL, MAY-JUNE 1969
frequency analysis, and time interpolation of data from one sampling
rate to any other sampling rate. These applications are explained in
detail. The chirp 2-transform algorithm also permits use of a radix 2
fast Fourier transform program or device to compute the discrete
Fourier transform of an arbitraiy number of samples. Examples were
presented illustrating how the chirp z-transform algorithm was used
in specific cases. It is anticipated that other applications will be found.
VII. ACKNOWLEDGMENT
The authors would like to acknowledge the help and assistance of
Miss Carol McGonegal in programming many of the examples shown
in this paper and Miss Debbie Hougak for valuable clerical help in
preparing the manuscript.
APPENDIX
Fast Fourier Transforms for Two Real L Point Sequences
The purpose of this appendix is to show how the fast Fourier trans-
forms of two real, symmetric L point sequences can be obtained using
one L/2 point fast Fourier transform.
Let x„ and y n be two real, symmetric L point sequences with corre-
sponding discrete Fourier transforms X k and Y k . By definition,
x n = x L - n n = 0,1,2, ■•• ,L- 1;
Vn = Vt-n
it is easily shown that X k and Y k are real, symmetric L point se-
quences, so that
X k = X h . k k = 0lt2t ... ,L- 1.
Y k = Y L - t
Define a complex, L/2 point sequence Un whose real and imaginary
parts are
Re [U„] = X 2n - ?/2n + l + V2n-1
Im [u„] = y 2n + x 2n+1 — x 2n -i.
n = 0, 1, ••• ,1/2-1.
The L/2 point discrete Fourier transform of u n is denoted U k and is
calculated by the fast Fourier transform. The values of X k and Y k may
CHIRP Z-TRANSFORM 1291
be computed from U k using the relations
X k = %{Re[U k ] + Re[U L/2 . k ]\
1
a 2ir ,
4 sin -=- k
Li
[Re[U k ] - Re[U L/2 . k }\
Y k = i{Tm[U k ] + lm[U W2 . k ]\
1
4 sin -=- k
Li
\lm[U k \ - lm[U L/2 . k )}
for k = 1,2, ••■ ,L/2 - 1.
The remaining values of X k and Y k are obtained from the relations
n-0
L-l
Y =TiVn
n-0
x L/2 = S (- dx
n-0
REFERENCES
1. Cooley, J. W. and Tukey, J. W., "An Algoritiim for the Machine Calculation
of Complex Fourier Series," Mathematics of Computation, 19, No. 90
(April 1965), pp. 297-301.
2. "What is the Fast Fourier Transform?," G-AE Subcommittee on Measure-
ment Concepts, IEEE Tran. Audio and Electroacoustics, AU-16, No. 2,
(June 1967), pp. 45-55.
3. Stockham, T. G., "High Speed Convolution and Correlation," 1966 Spring
Joint Computer Conference, Amer. Federation of Inform. Processing Soc.
Proc, 28, Washington, D. C, April 1966, pp. 229-233.
4. Helms, H. D. "Fast Fourier Transform Method of Computing Difference
Equations and Simulating Filters," IEEE Trans. Audio and Electroacous-
tics, AU-16, No. 2 (June 1967) , pp. 85-90.
5. Bluestein, L. I., "A Linear Filtering Approach to the Computation of the
Discrete Fourier Transform," 1968 Northeast Electronics Research and
Engineering Meeting Record, 10, November 1968, pp. 218-219.
6. Rader, C. M. and Gold, B., "Digital Filter Design Techniques in the Fre-
quency Domain," Proc. IEEE, 55, No. 2, (February 1967), pp. 149-171.
1292 THE BELL SYSTEM TECHNICAL JOURNAL, MAY-JUNE 1969
7. Schafer, R. W. and Rabiner, L. R., Automatic Formant Analysis of Speech
Using the Chirp z-Transform Algorithm" to be presented at the 1969 IEEE
International Conference on Communications, Boulder, Colorado, June,
1969.
8. Gold, B. and Jordan, K. L., "A Direct Search Procedure for Designing Finite
Duration Impulse Response Filters," IEEE Trans, on Audio Electro-
acoustics, AU-17, No. 1 (March 1969) pp. 33-36.
9. Gentleman, W. M. and Sande, G., "Fast Fourier Transforms for Fun and
Profit," 1966 Fall Joint Computer Conference, American Federation of
Information Processing Societies Proc, 29, Washington, D. C, November
1966, pp. 563-578.